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The seismic processes are well known to be self-similar in both spatial and temporal behavior. 
At the same time, the Burridge-Knopoff (BK) model of earthquake fault dynamics, one of the basic 
models of theoretical seismicity, does not posses self-similarity. In this article an extension of BK 
model, which directly accounts for the self-similarity of earth crust elastic properties by introducing 
nonlinear terms for inter-block springs of BK model, is presented. The phase space analysis of the 
model have shown it to behave like a system of coupled randomly kicked oscillators. The nonlinear 
stiffness terms cause the synchronization of collective motion and produce stronger seismic events. 
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I. INTRODUCTION 

The earthquake - a sudden stress relief in earth crust - is believed to take place when and where the stress exceeds 
certain critical value. At the same time, it is well known that the rock materials are not continuous in physical sense: 
they consist of grains or crystals of different sizes separated from each other by cracks of all possible sizes. The 
picture is more or less the same from the typical grain sizes up to the external size of the seismic zone. The situation 
very much resembles that in developed turbulence, described by Kolmogorov [Kol41a]. Therefore, we have to deal 
with a system with well manifested self-similarity properties. The self-similarity of seismic processes has got a lot 
of attention in phenomenological studies (see e. g. [Kag94,KB94] and references therein). A theoretical account for 
self-similarity of earth crust fracturing was given by Newman and Knopoff jNK90| , |Kol4]H and some other authors 
on the base of renormalization approach. The dynamical modeling of self-similar velocity weakening friction in BK 
model was recently presented in [3VR96|. Their hierarchical BK model accounts for cascade rupture propagation and 
seems to give a simple explanation to the Coulomb friction law so widely observed in nature. However, their model 
does not account for nonlinear elasticity which is also significant for earth crust deformations [ ZL95[ . 

The crux of our paper is the direct account of self-similarity (and scale-dependence, as it will be shown below) in the 
spirit of self-similar elasticity [ZL95|. Taking self-similarity as a given property of the rock material, we have changed 
the harmonic springs (U = ^§-) in BK model flBK6?H to nonlinear ones possessing self-similarity (to be explained 
below). We compare our simulations with that of standard B K mode l CL89 , [SVR93|| and found our model to be more 
realistic in the sense of foreshokes and aftershocks clustering | Ger93 |. Besides that, the nonlinear elasticity, which we 
have yielded from the self-similarity, leads to the nonlinear stiffness of the form 



f(x) — —k(x)x, where k(x) = ko 



1 



where o is the length of relaxed spring. This nonlinearity is adequate to well known empirical fact of nonlinear stiffness 
of the crust. 

The remainder of this paper is organized as follows. In section 2 we remind the original BK model and analyze 
its virtues and shortcomings. The self-similar stiffness modification of BK model is presented in section 3. Section 
4 is devoted to the numerical algorithm used for simulations. In section 5 we compare our model and the original 
Burridge-Knopoff one. The statistical characteristics - autocorrelation function and the Hurst exponent - are also 
presented here. In section 6 we give some dynamical analysis of the model, which shows that the BK system behaves 
like a system of kicked oscillators. 



II. THE BURRIDGE-KNOPOFF MODEL 



There is no mechanical model capable of simulating all the features of real seismic process. However, to study 
the temporal characteristics of earthquake dynamics is sufficient to isolate a few main properties of the process. The 
feature isolated and put as a basis for the mechanical model of earthquakes was a conjecture that earthquake faults are 
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retarded by nonlinear friction between blocks of rock material. Such a purely dynamical model obeying the Newton 
laws, without any ad hoc taken random forces, has been proposed by R. Burridge and L. Knopoff 30 years ago [BK67|. 

The BK system is a system of N blocks of mass m^i — 1,N rested on a rough surface and connected by harmonic 
springs of stiffness k c to each other; each block is attached by a leaf spring of stiffness k p to moving upper line, see 
Fig. 1. 




FIG. 1. Fhe geometry of the BK model: the system is composed of N identical blocks of mass rrn 
"horizontal" springs, k p is the stiffness of pulling springs, v is the constant velocity of pulling line. 



k c is the stiffness of 



Initially (at t = 0) the system is at rest, and the elastic energy accumulated in "horizontal" springs is only due to 
randomly generated small initial displacements of the blocks from their neutral positions. The moving upper line, 
which simulates the movement of external driving plate, exerts the force f n = —k p (x n — vt) on each n-th block. The 
nonlinear friction is defined in such a way that it holds each block at rest until the sum of all forces applied to this 
block exceeds certain critical value Fq. Then the block makes a slip inhibited by nonlinear friction to a new position. 
A pause between two slips is believed to account for a pause between earthquakes. 

Despite its simplicity, the BK model enables to simulate sequences of slipping events, similar to that of real 
earthquakes. The distribution of event size, generated by BK model also resembles the Gutenberg- Richter power law 
]GR44 1 . However, for practical purposes of generating artificial earthquake catalogs, the "reliability" of BK toy model 
is not always sufficient. A number of attempts have been made to modify BK model to gain more realistic sequences 
of events, to gain more sharp clusterization first of all [ Ger93| . 

Following ]CL89| we rewrite the equations of motion in using dimensionless coordinates: 



= l 2 {Un+l - U n ) + l 2 (u n -i - U n ) - U n + VT + F(u n ), 



(1) 



where u n = ^f" is normalized displacement, I = yjk c /k p is normalized sound velocity, v = -j~p- is dimensionless driv- 
ing velocity, Fq is the amplitude of friction force. The upper dots denote differentiation with respect to dimensionless 
time t — ui p t, 



= The adimensional version of friction force [BK67] is taken in the form 



H 



F(u) 



F 



1+Dtgt-H) 



1-D(ii+H) 



for \ll\ < H 

- E(u - H) for u> H 

- E(u + H) for u<-H 



(2) 



where E is a linear friction coefficient, H is a threshold velocity below which blocks can move without causing fracture 
effects; I? is a parameter of friction nonlincarity. 

To keep with [BK67 CL89| we have chosen E = Fq = k p = m = 1,D = 5, v = 0.01. The velocity threshold value 
was taken H = 0.001. We have performed simulations for 10, 50 and 100 blocks chains. 

The BK system has several evident time scales. The characteristic time T p = uj p l is the typical time of harmonic 
oscillation; tl = 1/v (dimensionless units) is the time taken by driving force u — vt to reach the critical value Fq = 1 
and to cause a slip; t s = N/l is the time required by a "sound" wave traveling with velocity I to pass over entire chain 
of N blocks. The ratio of the two latter parameters 9 — is a n important parameter of the system, which accounts 
for the transition from chaotic to solitary wave regime |SVR93 SVR96|. However, in the discrete spring-block model 
(Q) there is no dependence on the relaxed spring length a. It seems natural on one hand - the seismic processes are 
scale-invariant practically at all scales, - but on the other hand the model admits the possibility of block overlapping, 
when simulated numerically. The latter is completely unphysical. 
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III. SCALE DEPENDENT BK MODEL 



A. Scale invariant stiffness model 

Regardless all the merits of spring-block model evidently it does not adequately describe the elastic properties of 
rock material. 

Let us consider a spring of relaxed length a. If we contract or dilate it by Ax, we exhibit a reaction force / = —kAx] 
where k is its stiffness. However, it is absolutely impossible to contract a spring up to zero length, Ax — a. More 
than that, as any geologist knows, in practice k is not the same for large contractions and dilations. More precisely, 
the stiffness monotonously decreases when going from contractions to dilations, and then, achieving certain critical 
value, rock crashes. 

The simplest way to describe such behavior is to add cubic and quartic terms to spring potential energy 

[/(Ax) = + A (Ax) 3 + B(Ax) 4 . 



However, since we are going to deal with earth quake processes, we also have to account for such a well manifested 
property of the crust as self-similarity. Following ||ZL95 l, let us take the stiffness to be homogeneous of degree e 



fe(Ai) = \- e k(l), (3) 
normalized to the rigidity of relaxed spring k(a) — kg. The energy of the finite deformation Ax is then given by 



U(Ax) = 

After straightforward calculation this gives 

U(u) = k a 2 — " — 



Ax r Ax 

k(y)ydy = 
o Jo 



a"k ydy 



-(1-(1- U ) 2 -) 



(4) 



(5) 



where u = £z is the relative deformation. The graph of the function U{u) is presented in Fig. 2. 
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FIG. 2. The graph of the potential energy of a spring with self-similar stiffness k(Xa) — X t k(a), function (jE]). X-axis values 
correspond to relative deformation u = — , where a is the length of relaxed spring. Plotted for e = 0.5, fco = a — 1. 
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We will proceed with a case of small nonlinearity, which is the restriction of general formula (g) for the case of small 
e, and has convenient linear stiffness k = const as a limiting case at e — > 0. We will concentrate ourselves on elastic 
nonlinearity without taking friction into account at this point. Taking the Taylor expansion of (||) with respect to e 



U(u) = k a* [- + e( y + - + -) + e 2 (- + -) + o(e ») 



(6) 



In the first order with respect to e we get f(x) = —W- = —k\[x)x, where 



is the deformation dependent spring strength. Therefore, the self-similarity of elastic media necessarily causes non- 
linearity. 

Using the nonlinear potential @ for small displacements (ti < 1) we have 



U{u) = k^a 2 



(7) 



as a spring potential, we account for small influence of self-similar stiffness (e -C 1). 

The first term of the equation (Q), the harmonic spring potential, is already contained in the BK model. The 
nonlinearity comes from the second term. The force applied to the n-th block of the spring-block chain due to this 
nonlinear term is 



d 



3a " 8a 2 n ' 



(8) 



where A„ = x n — x n+ i. So, the nonlinear stiffness term will appear at the r. h. s. of equation of motion (|l|): 

1TlX n — k c (x n +\ 2x n -f- X n — \) kp{x n vt) -\- fnlin(x n ) -\- F{x n ) : 



where f n iin(x n ), given by (§ 
above takes the form 



In dimensionless form the spring-block equation (Q) with nonlinear terms described 



I 2 (Un+l 



+ a[(u n 



+i - u n 



u n ) - l 2 (u.„ 
2 - (u„ 



U n -l)) - U n + VT 

u n . x f] 



(9) 



+ (i[{u n+ i - U n f - (u n - U n -l) 

where the dimensionless nonlinearity parameters are 

l 2 eF 



F(u„), 







PeFl 
2a 2 fc2 : 



all other notations are the same as for the equation (Q). 



IV. NUMERICAL IMPLEMENTATION 



In this section we consider the Burridge-Knopoff model consisting on N blocks. For numerical simulations we 
rewrite equation ([!]) separately for boundary (i — 1, N) and internal blocks (1 < i < N) motion in the following form: 



x\ = l 2 (x 2 - Xi) - (xi - vt) + x 2 , 

xn = 1 2 (xn-i ~ xn) ~ {%n ~ vt) + ® N (x N ,x N ^i,x N ), 

Xi = l 2 {x i _i-2Xi + Xi_i)-{Xi-ut) + ^i{x i -i,Xi,X i+ i,X i ), 

i = 2,...,N-l, 



(10) 



where functions comprise nonlinear terms of r. h. s. of (Q). Free boundary conditions were set for the first and the 
last blocks. 

The problem is reduced to the solution of the nonlinear equations of motion ( |To| ) with given initial conditions: 
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,N. 



(11) 



Let Pi(t) — ±i{t). We rewrite the system (JlOjJllJ) in the following form 



Xi 


= Pi(t), 


i = l,...,N, 




Pi 


= l\x 2 


- Xi) - (xi - vt) 


+ ^i(xi,x 2 ,Pi), 


PN 


= 1 2 (XN 


-l - x N ) - (x N - 


vt) + ^ N (x N ,x N -i,p N ), 


Pi 


= l 2 (x.^ 


i - 2xi + - 


(xi - vt) + Wi(xi-i,Xi,x i+ i,Pi), 




i = 


2,...,N-1, 





(12) 



with initial conditions: 



Xi(t )=u°, pi(t ) = v°, i = l, . 



, N. 



(13) 



This scheme has been used for numerical simulations of both linear BK model (e = 0) and that with nonlinear 
terms (e > 0). The nonlinear friction was given by (||). The initial displacements u, were randomly generated with 
the amplitude 0.0001, the initial velocities were set to zero, we choose I 2 — 100. Let St be the time step t; = to + I5t. 
To get an approximate solution of the Cauchy problem ( |l2] , |l5| ) we use the implicit numerical scheme: 



x\ - x]- 1 



oMpl+p'r 1 ), % = !,. 



(pi 



St 

Pi 



, N, 



z— 1 > 



= 0.5(l 2 (x l 2 - 4) - x[ + l 2 {x l 2 ~ l - x[~ l ) - 
+ ^ + ~-l>\ 



St 



J-l _i— 1 „i-l-\ 



— — 0.5(l 2 (x l N _ 1 x l N ) x l N + l 2 (x l N ^ 1 x^ 1 ) x N •*) 



(14) 



(,Pi 



P 



i — IN 



St 



+ vt + ^i(xl_\a 
keep with [BI 

velocity threshold value was taken H — 0.001 



0.5(Z 2 (^_i - 2x\ + x\ +1 ) -x\ + l 2 ix\r_\ - 2x l r x + x ^\) - x \^) 
' ' ' 1 x l r + \J^), i = 2,...,N-l, 



with initial conditions (|l|). To keep with [BK67 CL89| we have chosen E = Fq = k p = m = 1, D = 5, z/ = 0.01. The 



V. RESULTS 
A. General outlook 



The dynamical behavior of both linear (^) and nonlinear ( |ic| ) models was studied in terms of potential and kinetic 
energy. We calculated the time dependence of potential U and kinetic T energies: 



Uc(r) 


Ar_1 
= 2 E ( u *+i 

i=i 


(r)-^(r)) 2 , 




- -o- £ (Ml4 

d i=l 


iW^,(r)f 


U m (r) 


1 N 

i=l 




Upot 


= *7 C + Ef m , 




T(r) 


1 N 

i=l 





Z 2 e^\ 



(15) 



(16) 
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The dependence of the energies on dimensionless time r is shown in Figs. 3 and 4, for linear (e = 0) and nonlinear 
(e = 0.2) model, respectively. 
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FIG. 3. Potential and kinetic energies for the linear (e = 0) N=100 block BK system I = 10, F = 1, D = 5, H = 0.001, E = 1. 
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FIG. 4. Potential and kinetic energies for the nonlinear (e 
/ = 10, Fa = 1, D = 5, H = 0.001, E = 1. 



0.2), N=100 block BK system 



As it can be seen from the pictures, after the initial period of periodic slip-stick motion of the system as a whole 
(the typical period of one cycle is about loading time tl = tl = 100 for our simulations), the process becomes 
unstable and the system enters the chaotic regime. 

As it was shown in the original paper of R. Burridge and L. Knopoff, the spring block system described by the 
equation (|I]) has a trivial solution Vn : u n (t) = vt+F(y). This solution is unstable with respect to small perturbations. 
The time of instability of this uniform motion is comparable to the loading time tl , after which the system exerts a 
slip. 

In our numerical investigations we have found another type of instability. After a period of time the quasi-periodic 
motion, when the system slides and sticks as a whole, also becomes unstable and actually chaotic when most part of 
inter-block springs are exited starts. The typical time of this instability is approximately the same for both linear 
and nonlinear models, see Figs. 3 and 4. At this stage (r < 500) the inter-block springs k c have not accumulated 



7 



enough energy and the dynamic of the system is determined by nonlinear friction force. Later, when nonlinear terms 
have got sufficient energy and the role of cubic and quartic terms in ( |lfj| ) becomes significant, the nonlinear model 
shows stronger clusterization of events than the linear one. This fact was traced on the autocorrelation function (see 
Fig. 5.), 

_ (x(t)x(t + z)) 

V[Z) ~ (x^t)} • 1 ' 

The averaging was taking over the available time series of energies < r < 2000. The increasing of correlations at 
intermediate times (z < 3tx) is clearly observed for nonlinear model. 




20 40 60 80 100 120 140 



z 

Autocorrelation functions for the potential energy, Uc 



j i i 

20 40 60 80 100 120 140 

z 

Autocorrelation functions for the kinetic energy, T 



FIG. 5. Autocorrelation functions for potential and kinetic energies. Shown for the nonlinear e = 0.1, (dashed line) and linear 
models. N = 10 blocks for both curves 

This leads to a collective behavior and probably provides an inverse energy cascade - an energy flux from high 
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frequency modes to low frequency ones, like that observed, say, in hydrodynamic turbulence. Besides that we can say 
the system to mimic the liquification processes, taken place for real earthquakes. For strong earthquakes, long before 
the main earthquake event the foreshock motion of the crust are governed by elasticity equation; later - in co-seismic 
stage - we have more rheological behavior rather than elastic process. The gravity waves like that on shallow water 
can be observed on the ground. The origin of these phenomena is that after a period of foreshocks the crust becomes 
less connected, less continuous and behaves as a system of particles with Van-der-Vaals interaction |Lom94|. In this 
sense, the instability we have observed at r about 1000 (y = 0.01,7V = 10) is a solid-liquid phase transition, see 
Figs. 3,4. 



B. Clusterization 

The time dependence of potential and kinetic energies (H© for linear (e = 0) and non-linear (e = 0.2) model are 
presented on Figs. 3 and 4, respectively. 

The nonlinear potential terms @) contributes to "horizontal" spring energies (second graph on both pictures). This 
potential terms, being relatively small in absolute values, play an essential role in nonlinear effects, happening in 
the spring-block chain. Initially, up to the time t < 500, the systems (both linear and non-linear) make a periodic 
slip-stick behavior as a whole. The potential energy of horizontal springs is very low at this stage: the internal degrees 
of freedom are not yet activated. Approximately at r = 500 this periodic motion breaks down to quasi-periodic weakly 
chaotic regime. 

At this stage the nonlinear system shows stronger clusterization, than its linear counterpart; the groups of events 
looks more distinguished for nonlinear system, see Fig. 4. 

It can be seen from both potential energy of "horizontal" springs and the total potential energy (the graphs in the 
bottom). Possibly, for some other values of parameters a stronger intermittency will be observed. This, however, will 
be the subject of further investigation. Here we have just to note that nonlinear terms provide more coherence of 
seismic events (which can be observed as long wave modulation) and make the events stronger since more potential 
energy is accumulated. 



C. Solitons 



The dependence of displacement of coordinate and velocity as function of time and number of blocks for the case 
of nonlinear model e = 0.1, N = 10 is presented on Fig. 6. 
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coordinate 



FIG. 6. The displacements u (lower picture) and velocities u (upper picture) as functions of time and block number. Calcu- 
lated for nonlinear model of N = 10 blocks, e = 0.1 



The solitary waves moving from left to right are clearly seen on this pictures. 

The localized solutions (propagating kinks u(s,t) = u(r ± s/j3)) of the BK system were already mentioned by 
Carlson and Langer as solutions of continuous limit equation of motion (formula 3.6 of [CL89]): 



(1 - t 2 lf?)u = -u- <\){%v + 2Cu). 

In numerical simulations they were shown to coexist with chaotic modes. The properties o f solitar y wave solutions 
depend on sound velocity £ = al, dissipation parameters, and crucially (as it was shown in SVR93 | paper, specially 
devoted to solito ns in BK model) on the dimensionless momentum parameter O = Nv, later generalized by the same 
authors to the relation of loading time to traveling time 9. 

Roughly, the number of solitary waves was numerically found to be 0/8. In our model there is a special source for 
solitary waves: the nonlinear spring terms, arising from self-similarity. 
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To study the effects of nonlinearity in the continuous limit of BK system, let us rewrite (1) in the co- moving frame 
Vn = u n - vr: 



Vn = l 2 (Vn+l - 2y n + Vn-l) - Vn 



n+l 



Vn) 2 ~ (Vn ~ Vn-l) 2 



(18) 



+ (3[(y n +i - Vn) 3 - (y n - Vn-i) 3 ] + F(2(y n + 2(u). 



Following Zabuzsky and Cruskal [ZK65|, who studied solitons in a similar system of differential equations, let us 
suppose that at rest all the masses of relaxed system are located at points x n = na, labeled by block number n, then 
we can use the Taylor expansion for y n +i and y n -i- 



Vn±l 



Substituting this into (fL9|) we obtain 



y n ± ay n ' + -^y n 



i a 



4! ' 



± 



jjn = l 2 a 2 y n " - y n + 2aa 3 y n 'y n " + 3f3a A y n ' y n " + -^Vn"" + friction terms. 



(19) 



Up to the second and the forth terms at the r. h. s., the equation (19), we have driven at, coincides with the 



solitonic equation of the Cruskal and Zabuzsky paper [ZK65]. The latter, for the case of small a can be transformed 
to Korteveg de Vries like equation 



V = V 



2ey'y" 



3 ,2 „ 

r y y 



friction terms.. 



(20) 



The detailed analysis of this equation will be done in further investigations. 

However, we can stress either at this point that there is a principle difference between simple solitonic models of 
the paper [CL89| and our equation ([l9]). The former models come from no-scale approximation of BK model, which 
does not contain any characteristic scales. In this sense it much resembles the (Kolmogorov) universal regime of 
hydrodynamic turbulence pCo!41a . The equation ( |l9| ) comes from the finite grid approximation of an explicitly scale 
dependent (but self-similar) model, which is more adequate for seismic events with strong deformations. 



D. The distribution of events of different energies 



The original BK model was well tested to simulate empirical Gutenberg-Richter law [ GR44 ] , the relation between 
the probability of event and its seismic moment of the form 



log TV 



b log M, where b w 1 . 



In BK mode l simu lations the seismic momen t M was understood as a sum of displacements taken over all blocks 
M = J2i u i CL89 1 . In original paper |BK67] the linear dependence between the logarithm of cummulative number 
of events with seismic moment greater or equal to M and logM holds only for the events of intermediate magnitude. 
We have tested both linear and nonlinear model in the same way and get the same conclusion about the domail of 
its validity. The distribution of events of different magnitude is presented in Fig. 7. 
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FIG. 7. The logarithm of normalized number of events with seismic moment greater or equal to M vs. log M. Solid line is 
used for linear model, dashed line for nonlinear model of 10 blocks. 



E. Hurst exponents 



An important characteristics of any natural hazard process is the Hurst exponent of its strength [ Fed8? , Lom94 1 . 
That means the time power-law behavior of the maximal deviation versus dispersion ratio 



R(z) 
S(z) ' 

where R(z) is maximal deviation taken place for < t < z. 



(21) 



R(z) = max X(t, z) — min X(t,z), 

0<t<z 0<t<z 



where 



X(t,z)\ 



is the accumulated deviation for the same time period. In our model it will be accumulated potential energy. The 
square mean deviation for the process £(t) in the time domain < t < z is 



S(z) 



\ 



~E(£(m)-<£>*) 2 
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For Brownian motion, a purely random process, h = 1/2. For other stochastic processes, H can be less or greater 
than 1/2 (see [Fed88| for general explanation of this point). The Hurst exponents for kinetic and potential energy of 
N = 10 block (e = 0) system calculated up to z = 2000 are presented in Fig. 8. The calculations have been done in a 
straightforward way by the definitions given above. 
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FIG. 8. Hurst exponents for kinetic and potential energies. Linear system of 10 blocks 



The Hurst exponents h(z) may be found from the relations (£l|) for discrete time argument. For both linear and 
nonlinear models we can see from Fig. 8 that for the Burridge-Knopoff system H exponent is much greater than 
0.5. It lies about 0.6-0.8, stabilizing for lower values about 0.7, when time increases. This fact shows the persistence 
tendency of the simulating processes. This is very close to the value of Hurst exponent obtained for real earthquakes, 



which is also about 0.7 [Lom94] 
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VI. DYNAMIC DESCRIPTION OF BK SYSTEM 



To study the dynamical properties of BK model, i. e. the interaction of coupled blocks comprising the spring block 
chain, their effect on the motion of the system as a whole we rewrite the system (Q) in vector form: 



ii = -l 2 [A]u - (u - vte Q ) + F(u), 



(22) 



where [A] is [N x N] matrix 



1 -1 







[A] 



-12-10 








V o 










-1 



-1 



... 0^ 

• • • 

• • • 

-1 

2 -1 

-1 1 / 



and 

u= ujv) t , F = (F(ui) 

Let {Ai,ej} be the eigenvalues and eigenvectors of [A]: 

[A]e. 



,F(u N )f, e = (l,l,...,l) 5 



A,-e,: 



1,N, (ej,efe) 



(23) 



where Ai < A2 < . . . < A,, 



Using the representation u = Xh=i s i e i we rewrite the equations of motion of Burridge-Knopoff model (^) as 

Sj = -{X l l 2 + l)si + (eo,ei)vt+ (e,,F), i = 7^N, 



(24) 

where Si(t) — (u(f),ej). It should be noted that matrix [A] has exactly N different nonnegative eigenvalues Ai < 
A2 < . . . < Xn with Ai = 0. The eigenvector corresponding to the lowest eigenvalue is collinear to eo: 



ei 



1 



N 



eo- 



Thus, we reduce the original problem to the system of nonlinear oscillators interlinked by nonlinear friction force F. 
Using Green functions for one-dimensional Helmgholtz equation u + lo 2 u — and the orthogonality of basic vectors 
{e^}^ we express the solution of BK system (|22| ) in the form: 



N 

u(t) = vteo + e,; \a,i saxwit + hi cos turf + 



sin(u>i(t — x)) 



(ei,F(u(x)))dx 



(25) 



where u>\ = AJ 2 + 1. Constants and bi depend on the initial state of the BK model. 

As one can see from this equation, the behavior of the solution u(t) is governed by projections of dissipative force 
F onto the basis {e^}^ . In Figs. 9,10 the time dependence of projections (e,, F(u(x))), £ = 1,3 are shown for the case 
of linear model, e = 0. 
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FIG. 10. The projections of dissipative force on 3d and 4th eigenvectors of the system (fell). A 10 block system. 



The projection of the first eigenvector looks quasi-periodic, all others exhibit chaotic behavior, like a randomly 
kicked oscillator [Zas7C|. 

For the first harmonics the convolution in r. h. s. of (pq ) is difficult to evaluate, and we use another representation 
to study its behavior qualitatively. Namely, we use the Green function for the equation s = g(t). Doing so, we obtain 



u(t) = isteo + 5Zj_ 2 ei a>i sin Uit + bi cos u>it + J Q 



l(u!i {t — x)) 



(ei,F(u(x))dx 



+ e 1 (a 1 + bit + (ei,f Q (t - x)(F(u(x) - u(x) + vxe )dx) 



(26) 



The function G(t) = ((F(u) — u + i/teo),ei) is shown in Fig. 11. This function looks like the above considered 
projections (F,e2), (F, 63), etc. 
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FIG. 11. Time dependence of driving term G(t). 
The phase portrait (Si, Si) of the first harmonics si is presented in co-moving frame 

Si(t) = si(t) - ut(e ,ei) 

in Figs. 12, 13 at different times for linear (e = 0) and nonlinear (e = 0.1) models, respectively. At the absence of 
friction force all orbits are evidently that of harmonic oscillator, i. e. ellipses. It should be noted that the motion of 
first mode Si up to a constant multiplier coincides with the motion of center of masses in co-moving ut frame. 
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FIG. 12. The phase portrait (Si, Si) for linear system of 10 blocks. 
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FIG. 13. The phase portrait (Si, Si) for nonlinear system of 10 blocks. 

At sufficiently large times, above t > 3000, the nonlinear system departs away from quasi-periodic attracting set of 
the linear model (compare Fig. 12 and Fig. 13). The horse-shoe like cycles shown in Figs. 12,13,14 and 15 evidently 
have self-similar structure if zoomed at different scales. The next-to-first harmonics (S2) phase portrait of nonlinear 
model (Fig. 15) is, in contrast, compressed in comparison to its linear counterpart Fig. 14. 
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FIG. 14. The phase portrait (S'2,5'2) for linear system of 10 blocks. 



20 



J I I I I I I I I I I I I I I I L 

-0.005 0.QQ5 



"0.06 tfcj | | |_| | | | | | | | | | | | |_j |_ 

-0.02 -0.01 0.01 



T = 1 100 



T = 2000 



0.3 






0.2 






0.1 









— ! \ 




-0.1 






-0.2 






-0.3 


—I 1 I 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 



-0.1 



-0.05 
T = 4000 



0.05 




T = SQ00 



FIG. 15. The phase portrait (£2, £2) for nonlinear system of 10 blocks. 
The higher modes s&, k > 2, have the same behavior. 

Therefore, we infer that the nonlincarity causes the energy transfer from individual block oscillations to the large- 
scale collective motions. This very much resembles gravity waves observed in co-seismic stage of strong earthquakes 
[Lom94]. These waves, looking like shallow water gravity waves are often developed on soft crust, where the nonlin- 
earity is much higher than on rigid, i.e. rocky crust. 

The representation of the solution u(t) in the form u(t) — allows to describe separately the behavior 

of the BK model as a whole and the interaction of blocks. This approach may be useful in futher studies devoted to 
the BK model. 
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CONCLUSION 



The dynamical behavior of BK mechanical model of earthquake faults was investigated. In contrast to standard BK 
we incorporate nonlinear terms in inter-block springs potential to account for self-similarity, which is widely observed 
for earth crust elastic properties. We have shown that nonlincarity arising from self-similarity can be considered as 
an additional source of solitonic behavior such as gravity waves observed in co-seismic stage of strong earthquakes. 
We have observed numerically that quasi-periodic slip-stick motion of the spring-block chain as a whole after a period 
of about 10 2 x tl breaks to chaotic behavior. This is much like a liquification process which plays an important role 
for earthquake caused disasters on soft ground. The phase analysis of the model shows sinchronization of phases of 
different blocks, which causes strong coherent motions of the system as a whole - large seismic events. 

So, being very simplistic the Burridge-Knopoff model imitates chaotization, related to the decay of low frequency 
modes to a number of high frequency ones, and the inverse processes - the formation of low frequency modes from 
differences of high frequency harmonics - known as inverse energy cascade. 
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